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In the Solar system the planets' compositions vary with orbital distance, with 
rocky planets in close orbits and lower-density gas giants in wider orbits. The 
detection of close-in giant planets around other stars was the first clue that this 
pattern is not universal, and that planets' orbits can change substantially after 
their formation. Here we report another violation of the orbit-composition 
pattern: two planets orbiting the same star with orbital distances differing by 
only 10%, and densities differing by a factor of 8. One planet is likely a rocky 
'super-Earth', whereas the other is more akin to Neptune. These planets are 
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twenty times more closely spaced — and have a larger density contrast — than 
any adjacent pair of planets in the Solar system. 

The detection of the first 'hot Jupiter' around a Sun-like star (1) was surprising to researchers 
who based their expectations on the properties of the Solar system. Theorists soon developed 
models of planetary 'migration' to explain how the planets' orbits can shrink (2). This led 
to a broader recognition that the architecture of planetary systems can change substantially 
after their formation. The planetary system reported in this paper is another example of an 
'extreme' planetary system that will serve as a stimulus to theories of planet migration and 
orbital rearrangement. The system features two planets in neighboring orbits with substantially 
different compositions. 

This system was discovered from the miniature eclipses or "transits" that cause the host star 
to appear fainter when the planets pass in front of the star. The target star (Kepler Object of 
Interest 277; hereafter, Kepler-36; also KIC 11401755, 2MASS 19250004+4913545) is one of 
approximately 150,000 stars that is subject to nearly continuous photometric surveillance by 
the Kepler spacecraft (3-5). The Kepler data revealed the transits of two planets. The loss of 
light during transits by planet b is only 17% as large as that from planet c (Fig 1). Both transit 
sequences show substantial timing variations (Fig 2). Indeed the variations are large enough that 
the smaller planet was initially overlooked by the Kepler data reduction pipeline (6, 7), which 
makes the usual assumption of strictly periodic orbits. It was subsequently identified with a 
search algorithm allowing for variations between successive transits up to a specified fraction 
of the mean period (8). 

The candidates have average orbital periods of P b = 13.84 d and P c = 16.23 d, which have 
a ratio close to 6/7. The transit timing variations are approximately linear in time, with abrupt 
slope changes after planetary conjunctions 6 x P b rs 7 x P c m every 97 d). The variations are 
anti-correlated; when one body is early, the other is late, and vice versa. This is diagnostic of 
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gravitational interactions (9-11), leaving no doubt that the two transiting bodies orbit the same 
star (as opposed to other scenarios in which two independent eclipsing systems are spatially 
unresolved by the telescope). 

Even without detailed modeling it is possible to show that the transiting bodies must be 
planets, by imposing the requirement of orbital stability. For long-term stability against colli- 
sions or ejections, the sum of the orbiting bodies' masses must be less than 10~ 4 times the mass 
of the star (12-14). This criterion leads to an upper bound of 40 Earth masses, hence planetary 
masses for both bodies. 

Precise knowledge of the masses and radii of the planets requires precise knowledge of 
those same properties for the star. This information is typically gleaned from spectroscopic 
observations. In this case there is an additional source of information: the Kepler data reveal 
solar-like oscillations of the host star, due to excitation of sound waves by turbulence in the 
stellar atmosphere. The frequencies of the oscillations indicate that the star has a density (25 ± 
2)% that of the Sun (15, 16). Analysis of high-resolution spectra of this star, subject to this 
density constraint, yields precise values for the stellar effective temperature and metallicity. 
The star is slightly hotter and less metal-rich than the Sun. This information combined with 
additional asteroseismic constraints (16) gives precise measures of the stellar mass and radius 
(Table 1). Based on these parameters (16), Kepler-36 is a 'subgiant' star, 2-3 billion years older 
than the Sun. 

We used the seismic results as constraints on the parameters of a photometric-dynamical 
model fitted to the Kepler data (16). This model is based on the premise of a star and two 
planets interacting via Newtonian gravity with an appropriate loss of light whenever a planet is 
projected in front of the star (17). Given initial positions and masses, the positions of the three 
bodies at the observed times were computed numerically, and the loss of light was computed 
assuming a linear limb-darkening law (18). Only data within half a day of a transit were fitted, 



4 



after correcting for extraneous trends by removing a linear function of time (16). 

We adjusted the parameters of the model to fit to the data. The optimal solution agrees well 
with the data, in particular with the transit-timing variations of both planets (Fig 1 and Fig 2). 
We calculated the credible intervals for the model parameters (16) with a Differential Evolution 
Markov Chain Monte Carlo algorithm (19) (Table 1). 

The masses and radii have uncertainties less than 8% and 3% respectively. This unusually 
good precision owes to the combination of the constraint on the planetary mass ratio (0.55 ± 
0.01) and planets-to-star mass ratio [(3.51 ±0.20) x 10~ 5 ] obtained from the detection of transit- 
timing variations, coupled with the usual geometric transit parameters and the asteroseismic 
constraints (16). 

Knowledge of the planetary masses and radii enables the study of their possible compo- 
sitions (16). Planet b's dimensions are consistent with a rocky Earth-like composition, with 
approximately 30% of its mass in iron. Any volatile constituents — those having relatively low 
boiling points such as hydrogen, helium or water — must be a small fraction of b's volume. As 
a limiting case, if the planet is assumed to have a maximally iron-rich core (20), water may not 
constitute more than 15% of the total mass, and any H/He atmosphere contributes less than 1% 
of the total mass. In contrast, planet c must be volatile-rich; it is less dense than water, imply- 
ing a substantial H/He atmosphere. Taking the interior to be free of water, with an Earth-like 
composition, the H/He atmosphere would contain 9% of the total mass. Even if the planet were 
half water by mass, the H/He atmosphere would still be > 1% of the total mass (Fig 3). 

As for the orbits, the photometric-dynamical model shows (16) they are nearly circular 
(eccentricities < 0.04) and coplanar (mutual inclination < 2.5 degrees). The evidence for 
coplanarity comes from the lack of significant variation in the transit durations (16,21). The 
orbits are also closely spaced: at conjunctions the planets approach one another within 0.013 
AU. The angular size of c as viewed from b would be 2.5 times larger than the full Moon viewed 
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from Earth. We investigated whether these close encounters are consistent with dynamical 
stability by scrutinizing a random sample of allowed model parameters (16). Direct numerical 
integration (16) showed that more than 91% of this sample avoided disruptive encounters over 
~0.7 million years. From this surviving population, 100 parameter sets were drawn randomly 
and numerically integrated for 140 million years; none experienced disruptive encounters. 

Of all the multiplanet systems for which densities of the planets have been measured, 
Kepler-36 has the smallest fractional separation between any pair of adjacent orbits, and this 
pair also has one of the largest density contrasts (Fig S20). These factors present an interesting 
problem. One would naturally suppose that the two planets formed at widely separated loca- 
tions in the protoplanetary disk, with volatile-poor b inside the 'snow line' at a few AU, and 
volatile-rich c outside. It will be interesting to see whether the usual 'migration' mechanism 
that is invoked to alter planetary orbits — tidal interactions with the gaseous protoplanetary 
disk — could draw together two planets from such different regions of the disk. Or whether 
the compositions and densities of the planets could have changed with time, due for example to 
the preferential erosion of the smaller planet's atmosphere by stellar irradiation (16) (Fig S21). 
Perhaps a combination of these factors will ultimately explain this puzzling pair of planets. 
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Star 

Mass, M* (M ) 
Radius, i?* (i? Q ) 
Mean Density, (g cm -3 ) 
Stellar Effective Temperature, T cS (K) 
Planet b 



1.071 ±0.043 
1.626 ±0.019 
0.3508 ± 0.0056 
5911 ±66 



Time of Transit, T b (BJD) 
Period, P b (day) 

Orbital Semimajor Axis, a b (AU) 
Mass, M b (M e ) 
Radius, R b (R®) 
Mean Density, p b (g cm" 3 ) 
Equilibrium Temperature, T cq b (K) 
Planet c 



2454960.9753^-^ 

h0.00082 



0.00060 



13.83989^ 

0.1153 ±0.0015 

A a c;+ - 33 

1.486 ± 0.035 
7 4fi+ - 74 

978 ± 11 



Time of Transit, T c (BJD) 
Period, P c (day) 

Orbital Semimajor Axis, a c (AU) 

Mass, M c (M e ) 

Radius, i? c (it!©) 

Mean Density, p c (g cm -3 ) 

Equilibrium Temperature, T cq c (K) 



2454955.9132^-^ 
Ifi 9Qsc;c;+o.ooo38 

lO.ZOOOO_ .00054 

0.1283 ±0.0016 

8-os±8:2g 

3.679 ± 0.054 
89+ - 07 

u -° y -0.05 

928 ± 10 



Table 1 : Characteristics of the Kepler-36 system. The stellar parameters are based on the 
analysis of the optical spectrum and the asteroseismic oscillations observed in the Kepler 
data (16). Because the orbits are not strictly periodic, the periods and times of transit given here 
refer to instantaneous ('osculating') values evaluated at an arbitrary reference epoch (2,454,950 
Barycentric Julian Date, BJD). The parameter ranges are based on the medians and 68% con- 
fidence limits of the marginalized posteriors. The equilibrium temperature was calculated as- 
suming a Bond albedo of A B = 0.3, with T eq = T cS jRj2a(l - A B ) l ' A . 
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Figure 1 : Kepler-36 light curve. Top. — Relative flux over 877 days, after applying a median 
filter to remove long-term trends caused by astrophysical and instrumental effects. The brief 
0.05% dips represent transits of planet c (also marked with red bars on the time axis). The indi- 
vidual transits of planet b (blue bars) are not easily discerned in the time series. The alternating 
pattern of transits (blue, red, blue,...) is flipped every w 97 days as the planets pass each other 
in their orbits. Bottom. — Composite transit light curves and best-fitting models, obtained by 
subtracting from each time stamp the nearest mid-transit time (calculated from the best-fitting 
model), and the data residuals. The gray bar shows the typical scatter of the relative flux. The 
transit of planet c (right panel, model in red) results in a loss of light about 6 times greater than 
that due to planet b (left panel, model in blue). 
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Figure 2: Visualization of transit- timing perturbations. Colored pixels encode relative flux 
(increasing from blue to green), with light gray pixels representing unavailable data. Each row 
represents an individual transit light curve, from the earliest observed transit at the bottom to the 
most recent transit at the top. The horizontal axis is the time modulo the mean orbital period for 
that planet. Strictly periodic transits would produce a blue vertical band. The curved, riverine 
bands that are observed are indicative of transit-timing perturbations. The red ticks on each 
row indicate the start and end time of transit, according to the best-fitting dynamical model. 
The range of relative fluxes in each panel spans the minimum and maximum data values in the 
middle panels of Figured] 
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Figure 3: Mass-radius diagram for small planets. Constraints for Kepler-36b and c are shown 
as two-dimensional joint probability densities and confidence contours (68% and 95%). Other 
exoplanets are shown for comparison: blue - the planets in Kepler- 11 (22), pink - Kepler- 
18b (23), gray - Kepler-20 b and c (24), brown - GJ 1214b (25), violet - CoRoT-7b (26), 
green - Kepler-lOb (27), orange - 55 Cnc e (28). Solar System planets are plotted using the 
first letter of their names (excluding Mercury). The curves represent theoretical models for 
planets of a given composition. Dotted curves are models of terrestrial bodies [those lacking a 
significant gas envelope (29)], with "Ice/Rock" - 50% ice and rock (silicates) by mass, "Rock" 
- 100% rock, "Earth-like" - 33% iron, 67% rock, "Iron" - 100% iron. Dashed curves are for 
planets with Earth-like solid cores surrounded with H/He envelopes with 5% or 10% of the total 
mass. The dash-dotted curve is for an Earth-like core and a water layer, in equal parts mass, 
surrounded by H/He envelope with 1.6% of the total mass. 
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Supporting Online Material 



In this supplement, we provide additional details regarding the discovery of Kepler-36 
(KOI-277, KIC 11401755, 2MASS 19250004+4913545). This supplement is organized as 
follows. In £pQ we describe the spectroscopic reconnaissance ( §1.11) and asteroseismol- 
ogy ( §1.21) of Kepler-36 that constrained the stellar properties, summarized in §1.3.31 In 
[JH we describe how constraints on the planetary orbits and bulk properties were in- 
ferred; we provide both a qualitative discussion ( ^2.11) and a more detailed description 
of the photometric-dynamical model ( §2.21) . model parameters, and the posterior estima- 
tion technique. In $3l we derive constraints on the composition of the planets based on 
their bulk densities and thermal histories. We discuss the short-term evolution in $4] and 
long-term stability of Kepler-36 in $5j We compare and contrast Kepler-36 with other 
known planetary systems in §[61 Finally, in § 111 we describe how to obtain the data used in 
this analysis and describe attached files containing additional information relevant to this 
work (including a sample of our MCMC results). 

1 Stellar properties of Kepler-36 

1.1 Spectroscopy 

1.1.1 Observations 

Spectroscopic observations to determine the stellar characteristics of Kepler-36 were conducted 
independently at two observatories. Two spectra were acquired using the Tull Coude Spectro- 
graph on the 2.7 m the Harlan J. Smith Telescope at the McDonald Observatory Texas on 29 
March 2010 and 8 April 2010. The two spectra were shifted and co-added to provide higher 
SNR for the spectroscopic analysis. One spectrum was acquired using the HIRES spectrometer 
on the Keck I 10 m telescope on 1 1 April 201 1. 

1.1.2 Determination of photospheric stellar parameters 

We used the Stellar Parameter Classification (SPC) method (30), to derive the stellar atmosphere 
parameters from the observed spectra. SPC cross -correlates the observed spectrum against a 
grid of synthetic spectra drawn from a library calculated by John Laird using Kurucz mod- 
els (31). The synthetic spectra cover a window of 300 A centered near the gravity- sensitive Mg 
b features and has a spacing of 250 K in effective temperature, 0.5 dex in gravity, 0.5 dex in 
metallicity and a progressive spacing in rotational velocity starting at 1 km s _1 up to 20 km 
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s . To derive the precise stellar parameters between the grid points, the normalized cross- 
correlation peaks were fitted with a three dimensional polynomial as a function of effective 
temperature, surface gravity and metallicity. This procedure was carried out for different ro- 
tational velocities and the final stellar parameters were determined by a weighted mean of the 
values from the spectral orders covered by the library. 

Initial values for the effective temperature and metallicity were determined and used as 
input in the asteroseismic analysis detailed in Section [L2l to estimate the surface gravity of the 
host star. SPC was then re-run, fixing the surface gravity to the value from the asteroseismic 
analysis (log(g) = 4.045). The final stellar parameters reported here, effective temperature T eff = 
5911 ±66^, metallicity [m/H] = —0.20 ±0.06 and vsini = 4.9 ±1.0 km -1 , are the mean values 
derived from the co-added McDonald spectrum. The HIRES spectrum was provided late to this 
analysis, but its parameters agree with those derived from the McDonald spectrum. 

1.2 Asteroseismic analysis of Kepler-36 

The asteroseismic analysis of Kepler-36 was performed by the Kepler Asteroseismic Science 
Operations Centre (KASOC) team. 

1.2.1 Estimation of asteroseismic parameters 

The analysis was based upon 15 months of Kepler short-cadence (1 -minute sampled) data, col- 
lected between 2010 March 2 and 2011 June 26. Fig. |ST] plots the frequency-power spectrum 
of the light curve (the planetary transit signals having first been removed). It shows a clear 
pattern of peaks due to solar-like oscillations that are acoustic (pressure, or p) modes of high 
radial order, n. The observed power in the oscillations is modulated in frequency by an enve- 
lope that has an approximately Gaussian shape. The frequency of maximum oscillation power, 

— 1/2 

z/ max , has been shown to scale to good approximation as gT eS 1 (32, 33), where g is the sur- 
face gravity and T c q is the effective temperature of the star. The most obvious spacings in the 
spectrum are the large frequency separations, Az/, between consecutive overtones n of the same 

1 /2 

spherical angular degree, I. These large separations scale to very good approximation as (p) ' , 
(p) oc M/R 3 being the mean density of the star with mass M and surface radius R (34). 

We used four independent analysis codes to obtain estimates of the average large separation, 
(Az/), and z/ max , using automated analysis tools that have been developed, and extensively tested 
(35-38) for application to Kepler data (15). A final value of each parameter was selected by 
taking the individual estimate that lay closest to the average over all teams. The uncertainty on 
the final value was given by adding (in quadrature) the uncertainty on the chosen estimate and 
the standard deviation over all teams. Excellent agreement was found between the results. The 
final values for (Az/) and z/ max were 67.9 ±1.2 pHz and 1250 ± 44 /iHz, respectively. 
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Figure S 1 : Frequency-power spectrum of Kepler-36, showing a pattern of nearly equally-spaced 
overtones of solar- like oscillations. 

Use of individual frequencies increases the information content provided by the seismic 
data for making inference on the stellar properties. We therefore also applied "peak bagging" 
to the frequency-power spectrum to extract estimates of the mode frequencies. Two codes 
were applied, one using a pseudo-global approach with maximum likelihood estimation (39); 
and another which performed a global fitting using Markov chain Monte Carlo (MCMC) (40). 
As per the average seismic parameters, we found excellent agreement in the fitted frequencies 
between the two approaches. The frequencies, z/J s , and their associated uncertainties, cx°^ s , 
are listed in Table 11.2.11 They were taken from the MCMC analysis, which provides more 
conservative frequency uncertainties. 

1.3 Estimation of stellar properties 

The analysis was divided into two parts. 

1.3.1 Grid-based modeling with average asteroseismic parameters 

In the first part we used a grid-based approach, in which properties were determined by search- 
ing among a grid of stellar evolutionary models to get a best fit for the input parameters, which 
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Table SI: Estimated frequencies z/°|? s and uncertainties a°^ s (68 % credible region) of Kepler-36 
(in yuHz). 
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were (Av), v max , and the spectroscopically estimated T eS and [Fe/H] of the star. Descriptions 
of the grid-based pipelines used in the analysis may be found in {41-44). An important out- 
put of this first part was not only a guideline set of stellar properties - which could be used 
as the starting points for detailed modeling of the individual frequencies in the second part of 
the analysis - but also iterated, and improved estimates of T e s and [Fe/H]. Initial values of the 
spectroscopic data were used together with the average seismic properties to estimate log g. The 
spectroscopic analysis was then repeated with log g fixed at this value, to yield the revised (and 
final) values of T cS = 5911 ± 66 K and [Fe/H]= -0.20 ± 0.06. Convergence of the inferred 
properties (to within the estimated uncertainties) was achieved after this one iteration. 

1.3.2 Modeling with individual mode frequencies 

In the second and final part of the analysis, the individual frequencies z/°]? s were used as the 
seismic inputs to a detailed modeling performed by four members of the KASOC team (TM, 
JCD, AM and SB). The spectroscopic inputs were also used, as per the grid-pipeline analysis. 
Estimated stellar properties from the first, grid-based part were used either as starting guesses or 
as a guideline check for initial results. Each analysis sought to minimize a standard x 2 metric. 
A separate value of x 2 was calculated for the asteroseismic and spectroscopic constraints, and 
these values are averaged for the final quality metric. 

TM used the Asteroseismic Modeling Portal (AMP), a web-based tool tied to TeraGrid 
computing resources that uses a parallel genetic algorithm (45) to optimize, in an automated 
manner, the match to observational data (46,47).AMP employs the Aarhus stellar evolution code 
ASTEC (48) and adiabatic pulsation code ADIPLS (49). Models were made using the OPAL 
2005 equation of state and the most recent OPAL opacities supplemented by opacities at low 
temperature (50), nuclear reaction rates from (51), and helium diffusion and settling following 
(52). Convection was treated with standard mixing-length theory without overshooting (53). 

Each AMP model evaluation involves the computation of a stellar evolution track from 
the zero-age main sequence (ZAMS) through a mass-dependent number of internal time steps, 
terminating prior to the beginning of the red giant stage. Exploiting the fact that (Av) is a 
monotonically decreasing function of age (46), the asteroseismic age is optimized along each 
evolutionary track using a binary decision tree. The frequencies of the resulting model are 
then corrected for surface effects following the prescription of (54). The optimal model is 
then subjected to a local analysis that uses singular value decomposition (SVD) to quantify the 
uncertainties of the final model parameters (55). 

JCD applied a fitting technique that has been used for the analysis of the Hubble observa- 
tions of HD 17156 (56), and Kepler observations of the Kepler exoplanet host stars HAT-P- 
7 (35) and Kepler-10 (57). The stellar modeling was carried out with the ASTEC code (48). 
The calculations used the OPAL equation of state tables (58) and OPAL opacities at tempera- 
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tures above 10 4 K (59); at lower temperature the opacities (60) were used. Nuclear reactions 
were calculated using the NACRE parameters (61). Convection was treated using the (53) 
mixing-length formulation. Frequencies were computed for the models using ADIPLS (49). 
The prescription in (54) was again applied in an attempt to deal with the surface term. 

For each evolutionary sequence in the grid of ASTEC models, the model Ai' min whose 
surface-corrected frequencies provided the best x 2 match to the observations was selected. The 
best match was obtained from application of homology scaling, i.e., it was assumed that fre- 
quencies in the vicinity of M' min could be obtained as ru n i(M' min ), wherer = [R/ ' R(M' min )]^ 1 - 5 
R being the surface radius of the model. A best-fitting model was then determined by minimiz- 
ing the sum {^nf ~~ ru ni(-M' mhl )^ I ( a °d S ) over a ^ observed modes, as a function of r. 
The resulting value of r defines an estimate of the radius of the best-fitting model along the 
given sequence. The other properties of this best-fitting model are then determined by linear 
interpolation in R, to the radius of the best-fitting model. Statistical analysis of the ensemble of 
best- fitting properties from all evolutionary sequences then yielded the final stellar properties, 
and their uncertainties. A Monte Carlo simulation, involving 100 realizations of the observed 
frequencies with the addition of normally distributed errors having the inferred standard devi- 
ations, showed that the observational errors on the frequencies make a modest contribution to 
the errors in the inferred stellar properties. 

AM used stellar models computed with the CLES code (62), from the pre-main sequence up 
to the sub-giant branch. For each evolutionary track adiabatic oscillations were computed for 
about 120 main-sequence and sub-giant models with the LOSC code (63) Grids of models both 
without and with diffusion and settling of helium and heavy elements (see (62) were considered. 
The opacity tables are those of OPAL96 (59) complemented at T < 10 4 K with the opacities 
of (60). The metal mixture used in the opacity tables was the solar one, as per (64). The 
nuclear energy generation routines are based on the cross-sections by NACRE (61), with the 
OPAL equation of state (65). Convection transport was treated with the classical mixing length 
approximation of (53). 

The analysis was performed in a similar manner to JCD, including treatment of the surface 
term using the prescription (54), and interpolation in R close to the x 2 minimum to find the 
best- fitting properties. 

SB made use of the Yale stellar evolution code, YREC (66) to model the star. The input 
physics included the OPAL equation of state tables of (65), and OPAL high-temperature opac- 
ities (59) supplemented with low -temperature opacities from (60). All nuclear reaction rates 
were from (67), except for the rate of the 14 N(p, 7) 15 reaction, which was fixed at the value 
of (68). Models were constructed for two values of core overshoot, and 0.2H P . Two families 
of models were constructed, one that included the diffusion and settling of helium and heavy 
elements as per the the formulation of (69), and one that did not include any diffusion and 
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settling. 



YREC was used in an iterative mode. In this mode the final T c g and radius for a model of 
a given mass and metallicity was specified, and for a given mixing length parameter a the code 
iterated over the initial helium abundance Y until a model with the specified T cS and radius was 
found. Note that this is similar to the construction of standard solar models, though in that case 
the iteration is performed over both the mixing length parameter and Y (since the solar age is 
a known constraint). Since the age of Kepler-36 is not known independently, the iteration over 
Y was done for many different values of the mixing length parameter. All solutions for which 
the initial helium abundance was less than the primordial helium abundance, Y p were rejected. 
Y p was assumed to be 0.245. 

The surface term correction was handled in a manner that is slightly different from the other 
analyses. The first step was the construction of a standard solar model with exactly the same 
physics as that used to model Kepler-36. This yielded the set u n i Q of solar model frequen- 
cies. These were then used to estimate a set of "surface term" frequency offsets, 8v n i Q , for the 
Sun by computing differences between the solar model frequencies and the solar low-degree 
frequencies observed by the Birmingham Solar Oscillations Network (BiSON), as in (70). 

Then, for each stellar model M.' under consideration, v n \ & and 5v n i & were scaled to the 
mass and radius of M.' using the homology scaling r = (Av(AA')) / (Au n i :Q ). The resulting 
ri^niQ-rSuniQ relation was then used to correct the stellar model for the surface term. Using a 
least squares minimization a factor /3 was selected so as to minimize (vffi — z^° rr ) / (o"nl >s ) 




over all observed modes, where z/^° rr = v^' + f3 r5u n i & , with r5u n i Q evaluated at ru n i Q = i/°j? s . 
1.3.3 Results on stellar properties 

Kepler-36 turned out to be a non-trivial star to model. The star is quite evolved, and a fraction 
of the models that can potentially provide a reasonable match to the observations show mixed 
/ = 2 modes in the frequency range where estimated frequencies are provided, thereby giv- 
ing the potential to complicate the analysis (71, 72). Since none of the observed frequencies 
showed evidence of strong mixing, model frequencies identified as those with the highest iner- 
tia, which did not satisfy asymptotic behaviour, e.g. (73, 74), were eliminated from the analysis. 
Some modifications were also made to the (54) surface-term prescription. This prescription was 
developed for application to radial modes: However, since for Kepler-36 the observed radial- 
mode data in Table [L2TT1 does not extend to the highest frequencies, JCD adapted the method 
so it could be applied to the non-radial modes. The correction was based on an estimate of the 
large separation (Au) and an average frequency obtained by fitting to all the observed modes. 

In spite of the above potential complications, differences in the treatment of the surface term, 
and differences in the physics inputs, there was excellent agreement in the properties returned by 
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all four analyses. For each of M, R, log g, mean density (p), and age r, the standard deviation 
of the four estimates was significantly lower than the median formal uncertainty of the property. 

Our final properties are those of the modeler whose estimates lay closest to the median (and 
mean) over all modelers. For the uncertainties we took the chosen modeler's uncertainty for 
each property and added (in quadrature) the standard deviation of the property from the different 
modeling results. We find M = 1.071 ± 0.043 M , R = 1.626 ± 0.019 R , (p) = 0.3508 ± 
0.0056 gem" 3 , hgg = 4.045 ± 0.009 dex and r = 6.8 ± 1.0 Gyr. (We add that the initial, 
grid-modeling-based estimates agreed to within uncertainties) Kepler-36 has a luminosity 2.9 
times that of the Sun and an approximate distance of 470 pc. 



1.4 Photometric determination of the stellar rotation period? 

The full Kepler photometry demonstrates a significant modulation with a period of 17.20 ± 
0.2 days (Figure |S2|) . This period corresponds to the rotational period of Kepler-36 if this 
modulation is due to the "transit" of photospheric surface features. In fact, this measured period 
agrees well with the maximum rotational period; we find P ro t,max = 16.8 ± 4.3 days assuming 
the results from above of vsini* = 4.9 ± 1.0 km/s and R* = 1.626 ± 0.019. A 17.20 day 
rotational period is typical (75) for subgiant stars with the same B — V color as Kepler-36 
(estimated to be B — V ~ 0.55). We note that this period is close to the orbital period of Planet 
c (at 16.28 days). However, it is extremely unlikely that the star has been tidally synchronized 
with c given its mass and distance. The correspondence is likely coincidental. 



2 Constraints on Planetary Orbits, Masses and Radii 

In this section we describe in detail how we derive the properties of the planets and the orbits 
from the photometric data. 



2.1 Qualitative Considerations 

Mass ratios between the planets and the star may be quickly estimated based on the qualitative 
behavior of the variations of the times of transit relative to a constant period model [i.e., the 
functional form of the transit timing variations (TTVs); see Figure |S3l Tables |S2] and |S3), as 
described as follows. 

The planetary mass ratio can be estimated assuming conservation of energy. Since the 
planets interact gravitationally, if they do not experience a close approach, they each maintain 
nearly Keplerian orbits. This suggests that 

n „ (7.U..U,, <7.U,.U, 
^orbital ~ E b + E c k, « const. (1) 
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Figure S2: Periodogram of the Kepler photometry at low frequencies demonstrating a strong 
peak at 0.058 d _1 or 17.20 days. This may represent the rotational period of Kepler-36; it agrees 
with the maximum rotational period determined from v sin % and the radius of Kepler-36. 
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Epochs Ephemeris, T Period, P (days) 

(BJD-2,454,900) 



1 ? 


fin Q771 -i- n D7fin 


1 ^ 

io. 


orwo 1 n 
ou i o nz u. 


U t 


5, 6, 7, 8, 


130.1465 ±0.0576 


13. 


8174 ±0. 


0196 


10, 11, 12, 13, 15, 


199.2772 ± 0.0546 


13. 


8310 ±0. 


.0145 


16, 17, 18, 19, 20,21,22, 


282.3275 ±0.0169 


13. 


8648 ± 0. 


0044 


23, 24, 25, 26, 27, 28, 29, 


379.4499 ± 0.0368 


13. 


8673 ± 0. 


0090 


30,31,32, 33,35, 


476.5412 ±0.0107 


13. 


8187 ±0. 


0033 


37,38,39, 40,41,42, 


573.3318 ±0.0037 


13. 


8244 ± 0. 


0010 


44, 45, 46, 47, 48, 


670.1201 ±0.0046 


13. 


8423 ± 0. 


0016 


50,51,52, 53,54, 55,56, 


753.2429 ± 0.0044 


13. 


8645 ± 0. 


0010 


57, 58, 59, 60, 61, 62, 63, 


850.3652 ± 0.0044 


13. 


8512 ±0. 


0012 



Table S2: Linear ephemerides for Planet b, measured independently from the Kepler photome- 
try. The ephemeris listed is that for the first epoch in the set under "Epochs." 

This equation neglects the interaction energy of the two planets as well as the the planet masses 
relative to the star in computing the orbital periods. Since a b /a c (P b /P c ) 2 ^, 



M h M, 



c 



p 2/3 p 2/3 



+ const. (2) 



Perturbing this equation gives 



M < „ SP b (P^\ 

This equation assumes that the energy exchange between the planets is small compared to their 
orbital energies (i.e. that there are no close encounters). Investigating Figure 2, 5P b ~ 8.5 hr 
and 5P C « 5.5 hr, while P c /P b w 7/6 « 1.17, so M c /M b ps 2. 

The ratio of the sum of the planet masses to the star can be derived from the total amplitude 
of the TTVs and the period of circulation (76). The best-fit periods are P b = 13.85 days and 
P c = 16.23 days and have a proximity to resonance of e = 1 — (j + l)P b /(jP c ) w 0.005 where 
j = 6 in this case. 

This distance from resonance causes the conjunctions to drift in inertial space over N con j ps 
j~ 2 e _1 pa 5 conjunctions. As can be seen in Figure 2, it takes ps 5 conjunctions (e.g. from 150 
to 600 days when the TTV traces cross one another in the same direction) to complete a TTV 
cycle. 

The semi-amplitude of the sum of the TTVs is St ~ 7 hours. Using equation 29 from (7(5), 
this gives a total mass ratio of 



M b + M c Ste* in _ 5 

A* = U = \ ~FT~ ~ 5 X 10 . (4) 
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This equation assumes that period variations dominate the transit timing variations (which can 
be seen in the riverplot as the jumps at conjunctions are smaller than the TTV due to changes 
in period between conjunctions). It also assumes plane-parallel planets of much smaller mass 
than the star, and neglects some factors that are of order unity. 

Resonance is expected for e < j 1 / 3 /! 2 / 3 = 0.0024, which is not satisfied for the 6 : 7 
commensurability. The best-fit solution's (described in § 12.2.51) resonant argument does not 
librate. 

Table S3: Planet c mid-transit times, measured indepen- 
dently from the Kepler photometry. 



Epoch 


Mid-transit time 
(BJD - 2,454,900) 





55.9135 ±0.0042 


1 


72.1543 ± 0.0041 


2 


88.3951 ±0.0071 


3 


104.6409 ± 0.0034 


4 


120.8930 ± 0.0039 


5 


137.1427 ± 0.0030 


6 


153.3924 ±0.0028 


7 


169.6447 ± 0.0036 


8 


185.8876 ± 0.0040 


9 


202.1310 ±0.0043 


10 


218.3740 ± 0.0041 


11 


234.6149 ± 0.0059 


12 


250.8567 ±0.0058 


13 


267.1010 ±0.0043 


15 


299.5228 ±0.0051 


16 


315.7430 ± 0.0035 


18 


348.1813 ±0.0032 


19 


364.4032 ± 0.0027 


20 


380.5782 ± 0.0039 


21 


396.8023 ± 0.0041 


22 


413.0243 ± 0.0065 


23 


429.2448 ± 0.0046 


24 


445.4671 ± 0.0045 


25 


461.6911 ±0.0067 


26 


477.8979 ± 0.0056 


27 


494.1446 ±0.0048 


28 


510.3889 ±0.0053 


29 


526.6322 ± 0.0044 
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Table S3: continued. 



Epoch Mid-transit time 
(BJD - 2,454,900) 



30 


542. 


.8778 


± 





.0040 


32 


575. 


3683 


± 





.0039 


33 


591. 


,6204 


± 





.0030 


34 


607. 


.8700 


± 





.0056 


35 


624. 


,1196 


± 





.0033 


36 


640. 


.3718 


± 





.0035 


38 


672. 


,8535 


± 





.0038 


39 


689. 


.0915 


± 





.0029 


40 


705. 


.3274 


± 





.0041 


41 


721. 


,5641 


± 





.0037 


43 


754. 


,0018 


± 





.0029 


44 


770. 


,2221 


± 





.0028 


45 


786. 


,4410 


± 





.0039 


46 


802. 


,6580 


± 





.0031 


47 


818. 


,8763 


± 





.0033 


48 


835. 


0969 


± 





.0038 


49 


851. 


,2757 


± 





.0048 


50 


867. 


,5075 


± 





.0055 


51 


883. 


,7375 


± 





.0050 


52 


899. 


,9658 


± 





.0048 


53 


916. 


,1958 


± 





.0051 


54 


932. 


,4275 


± 





.0041 



2.2 Photometric-Dynamical Model 

We modeled the Kepler light curve of Kepler-36 using a dynamical model to predict the motions 
of the planets, and a transit model to predict the light curve. We performed an initial fit to the 
measured times of transit (given in Table s IS2l and IS3l . using the dynamical model alone, to seed 
a subsequent fit with the full photometric-dynamical model. 

2.2.1 Description of the model 

The "photometric-dynamical model" refers to the model that was used to fit the Kepler photom- 
etry. This model is equivalent to that described in the analyses of KOI-126 (77), Kepler-16 (77), 
Kepler-34 and Kepler-35 (78). 

The underlying model was a gravitational three-body integration. This integration utilized a 
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Figure S3: Transit times of the two planets as predicted from a fit assuming a linear ephemeris 
between each conjunction for Planet b and individual epoch transit times for Planet c. The error 
bars represent the uncertainty on each transit time. Due to the smaller transit depth for Planet 
b, the uncertainties are much larger than for planet c. The dotted vertical line indicates the time 
at which short cadence data began to be collected for Kepler-36; the short cadence data yield 
much more precise times of transit. The open diamonds are the best-fit times according to the 
photodynamical model (§ I2.2I) . The typical uncertainty in the predicted times are 13 minutes 
and 3 minute for the transits of planet b and planet c, respectively. 
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hierarchical (or Jacobian) coordinate system. In this system, r b is the position of Planet b rela- 
tive to the star, and r c is the position of Planet c relative to the center of mass of Planet b and the 
star. The computations are performed in a Cartesian system, although it is convenient to express 
r b and r c and their time derivatives in terms of osculating Keplerian orbital elements: instanta- 
neous period, eccentricity, argument of pericenter, inclination, longitude of the ascending node, 
and time of transit: P b c , e b c , i biC , uj b:C , T b c , respectively. We note that these parameters do 
not necessarily reflect observables in the light curve; the unique three-body effects make these 
parameters functions of time. The "time of transit," in particular, does not exactly correspond 
to an measured transit time; exact correspondence would only be seen if the orbit proceeded in 
a Keplerian fashion from the reference epoch. 

The accelerations of the three bodies are determined from Newton's equations of motion, 
which depend on r b , r c and the masses (79, 80). For the purpose of reporting the masses and 
radii in Solar units, we assumed GM Sun = 2.959122 x 1(T 4 AU 3 day~ 2 and R Sun = 0.00465116 
AU. We used a Bulirsch-Stoer algorithm (81) to integrate the coupled first-order differential 
equations for r 6 c and r bjC . 

We did not determine the spatial coordinates of all three bodies at each observed time. 
Instead, to speed computation, we recorded only the sky-plane projected separation between 
star and planet and the sky-plane projected speed of planet relative to star at the calculated time 
of transit (for a given epoch). The times of transit were solved for numerically as being those 
times when the projected separation between the star and planet was minimized. The result of 
these calculations was a collection of transit times t\ , impact parameters b\ and speeds v\ for 
each planet k e {b, c] and for epochs ik E Nk where Nk is the set of observed epoch numbers 
for planet k. The motion of the planet relative to the star is approximately linear in the sky-plane 
such that the projected separation as a function of time is, to good approximation, 

+ Ot) (5) 

for times near (a few transit durations) of the calculated mid-transit time. 

The approximate photometric model for the relative stellar flux, f(t), is then defined as 

/(*) = i-E E / A (^( t )' i?fe ' M )[ 1 + c Wr 1 -o.5 < * - < 0.5 (6) 

k n • v, I otherwise 

where X(z, r, u) is the overlap integral between a limb darkened star of radius i?* (such that 
the radial brightness profile is I(p/R*)/I(0) = 1 — u[l — \J\ — (p/R±) 2 ] with linear limb- 
darkening parameter u) whose center is separated by a distance z from a dark, opaque sphere of 
radius r. The function c(t) is a piecewise constant function giving the extra flux, relative to the 
stellar flux, that is included in the photometric aperture and is assumed to be constant in each 
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Kepler quarter. X(z, r, u) may be computed semi-analytically with available codes {18). This 
photometric model, assuming constant transit velocity, is faster to compute than calculating the 
positions at each photometric cadence and results in a negligible change in the quality of the 
model fit to the data compared to exact integration. This model does not include the "anoma- 
lous" brightening events that occur when the Planet c occults Planet b during a transit (82). No 
such events are predicted within the current observations. 

The continuous model f(t) is integrated over a 29.4 minutes interval centered on each long 
cadence sample. The continuous model is compared "as is" to the short cadence observations. 

2.2.2 Local detrending of Kepler data 

The Kepler light curve ("SAP_FLUX" from the standard fits product) for Kepler-36, spanning 
ten Quarters (877 days), is reduced to only those data within 0.5 day of any transit of either 
planet. Only long cadence data (29.4 minute cadence) are available for the first 406 days of the 
available Kepler data. Short cadence data (1 minute cadence) are available for the remaining 
456 days. Short cadence data are used when available. Transits are missing only when Kepler 
data are unavailable. Data are missing as a result of observation breaks during quarterly data 
transfers or spacecraft safe modes. 

Each continuous segment of data (being those data observed within 0.5 day) has a local 
linear correction divided into it. The parameters of this linear correction are found through an 
iterative process, as described as follows. In the first step, we masked the higher signal-to- 
noise transits of planet c and then performed a robust linear least-squares fit to each continuous 
segment. The data, having divided out this correction, were then fit with the photo-dynamical 
model with a nonlinear fitter (Levenberg-Marquardt). The best-fit model was then used to 
identify cadences in transit which were then masked and the fit and division was repeated. This 
process was repeated until the corrections converged to a sufficient tolerance. 

We also experimented with performing this correction inline with the posterior estimation 
routine (described in §2.2.61) but found negligible differences compared to leaving these fixed 
according to the converged linear correction found from the process described above. 

2.2.3 Specification of parameters 

The reference epoch was chosen to be t = 2, 454, 950 (BJD), near a particular transit (a some- 
what arbitrary choice). 

The model has 20 adjustable parameters. Two parameters are related to stellar constraints 
from asteroseismology and are subject to the priors discussed in §1.3.31 the stellar density 
times the gravitational constant, Gp*, and the stellar radius, R±. Two parameters are the mass 
ratios q + = (M b + M c )/M* and q p = M b /M c . Four parameters encode the eccentricities and 
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arguments of pericenter in a nonlinear way that reduces the complexity of the posterior topology 
(resulting in effectively linear correlations in these parameters): 



h- = 0.895 x eb cos out, — e c cos ou, 



c 



(7) 



h 



+ = e h cos Ub + e c cos u c 



(8) 
(9) 
(10) 



k 



= eb sm cob — e c sin ui c 



k + = eb sin ojb + e c sin u c 



In hindsight, we should have used a different linear combination of parameters, namely 
a b e b ± a c e c , where e ; = ej(cosa;j, sintUj). This is due to the fact that the transit timing con- 
strains most strongly the distance at closest approach and orbital phase of the planets at closest 
approach. This is a function of the epicyclic amplitude and phase which is proportional to c^ei 
for each planet. The differences are well constrained since these determine the closest approach 
for the planets, while the sums are poorly constrained. The only gain in re-executing our anal- 
ysis with this parameterization would be more rapid convergence to the target distribution (via 
the method described in §2.2.61) . Having already reached sufficient convergence, we did not 
pursue this. 

The remaining osculating parameters, 7 in total, are the periods P b , P c , the orbital in- 
clinations i c , the times of transit T b , T c and the difference between the nodal longitudes 
AO = O c — Cl b . The absolute nodal angle relative to North cannot be determined. 

Two more parameters are the relative radii of the planets: r& = Rb/R-k and r c = i? c /.R*. 
One parameter, u, parameterizes the linear limb darkening law for the star (described above). 
The final two parameters describe the width of the probability distribution for the photometric 
noise of the long cadence and short cadence observations, assumed to be stationary, white and 
Gaussian-distributed (ctlc and use - see the next section). 

Additionally, we may specify 9 more parameters describing the function c(t), introduced 
above, describing the relative extra flux summed in the aperture. The nine parameters specify 
the constant extra flux in each Kepler quarter. There are 10 quarters in total, but the absolute 
level of this contamination cannot be determined from the photometry alone, but, the levels rel- 
ative to a single quarter may be specified. We chose this reference quarter to be quarter 7 having 
observed in preliminary fits that this quarter's contamination fraction was near the median of 
the sample of 10 quarters. The addition of these degrees of freedom did not significantly change 
the distribution of the 20 primary parameters. 

We note that a discrete degeneracy exists between the relative nodal angle, Afi, and incli- 
nations, %b and i c . The light curve model is invariant under the transformation (AO, i c ) — > 
(—AO, 7r — ib, 7i — i c ). Consequently, we expect the marginalized posterior distribution of AO 
to be symmetric about AO = 0, and the distributions of ib, i c to be symmetric about tt/2. 
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2.2.4 Priors and Likelihood 



We assumed uniform priors in 16 of the 20 primary parameters described in the previous sec- 
tion excluding h+- and For these latter four parameters, we enforced uniform priors in 
eccentricities and arguments of pericenter. For these priors, the probability density obeys 



p(h+-,k+-)dh + -k + - oc p(e bjC ,u bjl 



e&e, 



de h ,e c du) b) uj c OC d eb ^ e< ,d u , b} uj c (11) 



e b e c 



The likelihood £ of a given set of parameters was taken to be the product of likelihoods 
based on the photometric data, assumed-Gaussian asteroseismology priors, and the weighting 
to enforce priors (described above): 



L oc <j LC exp 



2^c 



x a S c sc exp 



g (Aif c ) 2 



(12) 



x exp 



1 {AGp* 



2 V a G P * 





"l (AR*\ 2 ~ 


1 


x exp 


2 V *R* J 


x 




e b e c 



where AF i ' is the zth photometric data residual (either long cadence or short cadence), 
AGp*/acp+ and AR^/gr^ are the deviates between the asteroseismic constraints in density 
and radius, and o^csc are the width parameters describing the photometric noise for either 
long or short cadence data. These noise parameters have values that are within a fraction of 
a percent of the root-mean- square normalized flux as calculated using the Kepler light curve, 
excluding the data near any transit. 



2.2.5 Best-fit model 

We determined a best-fit model by maximizing the likelihood C as defined above. The maxi- 
mum likelihood solution was found by determining the highest likelihood in a large draw from 
the posterior as calculated with a Markov Chain Monte Carlo simulation as described in §2.2.61 
Figure IS4l shows the data folded on the best-fit period for both planets, as in Figure 2, as 
well as the analogous map for the best-fit photometric-dynamical model. The transit direction 
has been scaled by the ratio 7 : 6 such that the transits of the planets adjacent to one another 
correspond to nearly coincident times. The gravitational kicks imparted by the planets on one 
another are apparent after each conjunction. The purple dots mark the best-fit time of tran- 
sit at each epoch; the best-fit times are listed in Tables |S4] and [S5] The best-fit model has 
X 2 = 99640.7 for 99718 degrees of freedom. If we ignore the transits of b in the model (or 
equivalently we let the radius of b be zero) we find A% 2 = 498; this corresponds to a detection 
of b at 22 g. 
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Table S4: Planet b mid-transit times from the best-fitting 
photometric-dynamical model. 



Epoch Mid-transit time 
(BJD - 2,454,900) 
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434.9066 


28 


448.7681 


29 


462.6286 


30 


476.5201 


31 


490.3483 


32 


504.1788 


33 


518.0120 


35 


545.6733 


37 


573.3377 


38 


587.1575 


39 


600.9805 


40 


614.8053 


41 


628.6283 
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Table S4: continued.. 



Epoch 


Mid-transit time 




(BJD - 2,454,900) 


42 


642.4481 


44 


670.1250 


45 


683.9643 


46 


697.8070 


47 


711.6506 


48 


725.4916 


50 


753.2334 


51 


767.0986 


52 


780.9644 


53 


794.8334 


54 


808.7028 


55 


822.5694 


56 


836.4339 


57 


850.3660 


58 


864.2151 


59 


878.0656 


60 


891.9193 


61 


905.7728 


62 


919.6233 


63 


933.4727 



Table S5: Planet c mid-transit times from the best- fitting 
photometric-dynamical model. 



Epoch 


Mid-transit time 
(BJD - 2,454,900) 





55.9117 


1 


72.1520 


2 


88.3923 


3 


104.6379 


4 


120.8902 


5 


137.1401 


6 


153.3900 


7 


169.6425 


8 


185.8856 


9 


202.1295 


10 


218.3731 
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Table S5: continued.. 



Epoch Mid-transit time 
(BJD - 2,454,900) 



11 


234.6144 


12 


250 8568 

.7 V / • v7 / V7 v7 


13 


267 1016 


15 


299.5242 


16 


315.7447 


18 


348.1836 


19 


364 4059 


20 


380 5804 

-.7 v - V 7 • \j v7 V 7 A 


21 


396 8042 

-.7 '7 V7 • v 7 V 7 X^J 


22 


413.0258 


23 


429.2460 


24 


445 4680 

J- J- .7 • -L V 7 v7 V 7 


25 


461 6917 


26 


477 8975 

A II* ^7 1_7 1 v7 


27 


494.1442 


28 


510 3884 

/ -L V 7 • -.7 v - v 7 -L 


29 


526 6316 


30 


542 8771 

v^7 -L * 1 1 _L 


32 


575 3677 

v^7 1 \j * ' J V7 1 1 


33 


591 6200 


34 


607.8698 


35 


624.1196 


36 


640 3720 


38 


672 8539 

V7 1 • v7 / - 7 ' 7 


39 


689 0919 

V7 v..- ' 7 • V / '.7 -L ' 7 


40 

t -r\J 


70^ ^977 

1 \JO. OZ, 1 1 


41 


721.5644 


43 


754.0017 


44 


770.2221 


45 


786.4409 


46 


802.6579 


47 


818.8762 


48 


835.0968 


49 


851.2751 


50 


867.5069 


51 


883.7367 


52 


899.9649 
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Table S5: continued... 

Epoch Mid-transit time 
(BJD - 2,454,900) 

53 916.1948 

54 932.4265 



Figure |S5] shows a plot of the normalized data versus the best-fit photometric-dynamical 
model. The correlation is very nearly linear showing that the model is an excellent fit to the 
normalized data. 

Finally, we have created a set of panels of every transit with the best-fit model in Figures 
IS6l and [S7] Figure |S6] shows the long cadence data while Figure |S7] shows the short cadence 
data. We also have created a movie (Movie SI), attached with this supplement, which shows 
the transits of planet c as a function of time. 

The parameters of the best-fit model are listed in Table |S6] and derived parameters from 
that best-fit solution are listed in Table |S7| In §4.11 we describe the short-term evolution of the 
best-fit solution. 

2.2.6 Parameter estimation methodology 

We explored the parameter space and estimated the posterior parameter distribution with a Dif- 
ferential Evolution Markov Chain Monte Carlo (DE-MCMC) algorithm (19). In this algorithm, 
a large population of independent Markov chains are calculated in parallel. As in a traditional 
MCMC, links are added to each chain in the population by proposing parameter jumps, and then 
accepting or denying a jump from the current state according to the Metropolis-Hastings crite- 
rion, using the likelihood function given in §2.2.41 What is different from a traditional MCMC 
is the manner in which jump sizes and directions are chosen for the proposals. A population 
member's individual parameter jump vector at step i + 1 is calculated by selecting two randomly 
chosen population members (not including itself), and then forming the difference vector be- 
tween their parameter states at step i and scaling by a factor T. This is the Differential Evolution 
component of the algorithm. The factor T is adjusted such that the fraction of accepted jumps, 
averaged over the whole population, is approximately 25%. 

We generated a population of 60 chains and evolved through approximately 500,000 gen- 
erations. The initial parameter states of the 60 chains were randomly selected from an over- 
dispersed region in parameter space bounding the final posterior distribution. The first 10% 
of the links in each individual Markov chain were clipped, and the resulting chains were con- 
catenated to form a single Markov chain, after having confirmed that each chain had converged 
according to the standard criteria including the Gelman-Rubin convergence statistics and the 
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-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 
Folded time (days) Folded time (days) 



Figure S4: Plot of the transits of planet b (left) and planet c (right). Top is data; middle is best-fit 
photometric-dynamical model; bottom is the residuals after subtracting the model. Refer to the 
caption of Figure 1 for additional details. The greyscale ranges from 0.9994 (black) to 1.0001 
(white) on a common scale between the top four panels. 
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1.0004 



1.0002 - 



>< 1.0000 - 



0.9992 




Totality Transit 
Planet b 



0.9995 0.9996 0.9997 0.9998 0.9999 1.0000 

Model Flux 



Figure S5: Correlation between the photometric-dynamical model flux and normalized Kepler 
photometry (i.e. after removal of a local linear fit about the transit). The red points are data, 
sorted in flux, averaged in bins with an equal number of data samples. The blue line denotes 
exact correspondence between model and data. 
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t c = 55.9 



1.0002 
1.0000 
0.9998 
0.9996 
0.9994 

1.0002 
1.0000 
0.9998 
0.9996 
0.9994 

1.0002 
1.0000 
0.9998 
0.9996 
0.9994 

P 1.0002 

4: 1.0000 

> 0.9998 

ns 0.9996 

a> 0.9994 
DC 

1.0002 
1.0000 
0.9998 
0.9996 
0.9994 



0.9998 
0.9996 
0.9994 

1.0002 
1.0000 
0.9998 
0.9996 
0.9994 




t c = 130.1 



t c = 137.2 t c = 144.0 t c = 153.3 t c = 157.8 t c = 169.7 t c = 171.7 






t c = 185.7 



t c = 199.4 t c = 202.2 t c = 213.1 t c = 218.5 t c = 227.1 t c = 234.5 





t c = 240.9 



t c = 250.9 t c = 267.0 t c = 268.6 t c = 282.2 t c = 296.3 t c = 299.5 




IT 




t c = 310.2 



t c = 315.7 t c = 324.0 t c = 337.9 t c = 348.2 t c = 351.7 t c = 365.0 





■ •■ #•«*. 




t c = 380.0 



t c = 393.3 t c = 396.7 t c = 407.1 t c = 413.1 t c = 421.0 t c = 429.3 






t c = 434.8 



t c = 445.5 t c = 448.6 t c = 462.1 





-10 10 



-10 10 -10 10 -10 10 

Time Relative to Segment Center (hr) 



Figure S6: Long cadence Kepler data used in this analysis. Each panel shows data for a con- 
tiguous set of cadences near transits of either Planet b or Planet c. The time in listed above each 
panel is near the center time of each panel less 2,454,900 (BJD). The red curve is the best-fitting 
photometric-dynamical model. The time and flux axis is the same in all panels. 
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L = 478.0 L = 490.2 L = 494.2 L = 504.0 L = 510.5 



1 .0005 
1 .0000 
0.9995 
0.9990 

1 .0005 
1 .0000 
0.9995 
0.9990 



t c = 51 7.9 t c = 526.7 




1 .0005 
1 .0000 
0.9995 
0.9990 

1 .0005 
1 .0000 
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-10 10 -10 10 



-10 10 -10 10 

Time Relative to Segment Center (hr) 



Figure S7: Short cadence Kepler data used in this analysis. Each panel shows data for a con- 
tiguous set of cadences near transits of either Planet b or Planet c. The time in listed above each 
panel is near the center time of each panel less 2,454,900 (BJD). The red curve is the best-fitting 
photometric-dynamical model. The time and flux axis is the same in all panels (but larger in 
flux than Figure |S6] in order to account for the larger relative error in the short cadence data). 
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observation of a long effective chain length in each parameter (as determined from the chain 
autocorrelation). 

2.2.7 Characteristics of the parameter posterior 

The parameter values and derived values reported in Tables |S6l [SVJ beside the best-fit values 
(see §2.2.51) . were found by computing the 15.8%, 50%, 84.2% levels of the cumulative distri- 
bution of the marginalized posterior for each parameter. Figure |S8] shows two-parameter joint 
distributions between all parameters. This figure is meant to highlight the qualitative features of 
the posterior as opposed to providing quantitative ranges. The numbers in that figure correspond 
to the model parameters in Table IS6l with the same number listed as "Index." Some non-linear 
degeneracies still remain in the form of "bananas" shown in the joint distribution plane, but 
these appear to have been properly traversed and mixed according to the convergence criteria 
listed above. 

Figure |S9] shows the posterior distribution in the eccentricity and argument of pericenter 
planes (this figure is discussed further in A number of arguments are permitted, however, 
only small eccentricities are allowed. At 95% confidence, e b < 0.039 and e c < 0.033. 

The three-dimensional "mutual" inclination between the planets' orbits, which may be de- 
fined as an osculating term at the reference epoch as 

cos/ = sin 4 sin i c cos AQ + cos 2 C cos 4, (13) 

is constrained to within a couple of degrees. Figure ISlOl shows the distribution of /; Table IS7l 
gives confidence intervals in the mutual inclination. The constraint on the mutual inclination 
derives from the lack of significant transit durations variations in the transits of planet c. To 
show this, we examined the long term variation of the sky-plane inclinations of each planet 
based on a draw of 100 initial conditions from our posterior (short term evolution of the posterior 
is described further in § H]). The inclinations show circulation of a forced inclination about 
a free inclination. We isolated and fitted the secular component (over the time scale of our 
observations) so as to compute the transit duration variation due to the secular change for planet 
c as a function of mutual inclination. We plot this for our 100 initial conditions in Figure ISTTl 
This figure shows a clear correlation between transit duration change and mutual inclination. 
For a mutual inclination of ~2.5 degrees - at the tail of our marginalized posterior distribution 
in that parameter (see Figure IS 101 ) - the transit duration changes by about 5 minutes over the 
roughly 900 days of observation or approximately 0.1 minutes/epoch. We compared this rate 
with what can be constrained from the data: A simple transit model fit to individual transit 
epochs of planet c shows that the measured uncertainty in duration is approximately 10 minutes 
at each epoch. Fitting a linear model to the durations of 50 transits of c (the number that was 
observed), assuming zero transit duration variation, would yield a 1-sigma uncertainty in the 
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Index 


Parameter Name 


Best-fit 


50% 


15.8% 


84.2% 




Mass parameters 













Mean Density, (g cm -3 ) 


0.3531 


0.3511 


-0.0052 


+0.0053 


1 


Mass sum ratio, q + ( x 10 5 ) 


3.40 


3.51 


-0.14 


+0.24 


2 


Planetary mass ratio, g p 


0.553 


0.550 


-0.009 


+0.010 




Eccentricity parameters 










4 


Cosine Difference, h- 


0.00397 


0.00393 


-0.00027 


+0.00028 


6 


Cosine Sum, h+ 


0.004 


0.004 


-0.018 


+0.020 


10 


Sine Difference, k- 


0.0270 


0.0258 


-0.0018 


+0.0015 


12 


Sine Sum, k + 


0.027 


0.017 


-0.030 


+0.012 




Planet b Orbit 










9 


Orbital Period, P h (day) 


13.83988 


13.83989 


-0.00060 


+0.00082 


11 


Time of Transit, t b (days since t ) 


10.97981 


10.97529 


-0.00583 


+0.00548 


13 


Orbital Inclination, if, (deg) 


89.52 


90.01 


-0.71 


+0.69 




Planet c Orbit 










3 


Orbital Period, P c (day) 


16.23853 


16.23855 


-0.00054 


+0.00038 


5 


Time of Transit, t c (days since t ) 


5.91174 


5.91315 


-0.00097 


+0.00109 


7 


Orbital Inclination, i c (deg) 


89.76 


89.98 


-0.53 


+0.54 


8 


Relative Nodal Longitude, Af2 (deg) 


0.25 


0.06 


-1.54 


+ 1.48 




Radius Parameters 










15 


Stellar Radius, (i? ) 


1.619 


1.626 


-0.019 


+0.019 


16 


b Radius Ratio, Rb/R* 


0.00857 


0.00838 


-0.00017 


+0.00016 


17 


c Radius Ratio, R c / R* 


0.02079 


0.02076 


-0.00018 


+0.00018 


14 


Linear Limb Darkening Parameter, u 


0.446 


0.456 


-0.024 


+0.024 




Relative Contamination, F cont /F± (xlOO) 












Quarter 1 


1.41 


1.50 


-3.10 


+3.27 




Quarter 2 


-0.20 


-0.33 


-2.46 


+2.46 




Quarter 3 


1.75 


0.77 


-2.53 


+2.61 




Quarter 4 


-0.98 


-0.22 


-2.81 


+2.79 




Quarter 5 


2.34 


1.97 


-2.56 


+2.64 




Quarter 6 


3.43 


2.97 


-2.32 


+2.39 




Quarter 7 






(fixed) 






Quarter 8 


-0.02 


-1.41 


-2.31 


+2.37 




Quarter 9 


3.71 


3.76 


-2.27 


+2.28 




Quarter 10 


-0.72 


-2.29 


-2.02 


+2.12 




Noise Parameters 












Long Cadence Relative Width, a LC (xlO 5 ) 


7.89 


7.87 


-0.10 


+0.11 




Short Cadence Relative Width, a S c (x 10 5 ) 


33.692 


33.681 


-0.077 


+0.078 



Table S6: Model parameters as defined in the text. The reference epoch is t =2,454,950 
(BJD). 
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n non7o 
U.UzU / y 


U.UzU ( 


n nnm a 
— U.UUUlo 


1 n nnm a 
+U.UUUI0 


UllUdCL x dldlllCLCI Ul x IdllCL U, / JT^ 


n 1 ^ 


n nn 

— u.uu 


n 1 8 
— u. 10 




Impact Pdrameter of Pldnet c, b c / i?* 


0.07 


0.00 


-0.16 


+0.16 


Trdnsit Velocity of Pldnet b, v b /R* (ddy -1 ) 


7.129 


7.079 


-0.104 


+0.043 


Transit Velocity of Pldnet c, v c / (ddy -1 ) 


6.579 


6.541 


-0.087 


+0.035 


Transit Duration of Pldnet b (hr) 


6.74 


6.78 


-0.11 


+0.09 


Transit Duration of Pldnet c (hr) 


7.430 


7.443 


-0.020 


+0.020 


Transit Ingress/Egress Duration of Pldnet b (min) 


3.490 


3.455 


-0.086 


+0.122 


Transit Ingress/Egress Duration of Pldnet c (min) 


9.12 


9.21 


-0.14 


+0.27 


Temperature Scdling of Pldnet b, \J R+/2a b 


0.18087 


0.18104 


-0.00045 


+0.00045 


Temperature Scdling of Pldnet c, R±/2a c 


0.17149 


0.17165 


-0.00043 


+0.00043 



Tdble S7: Derived Pdrameters. The reference epoch is t =2,454,950 (BJD). 



linear model slope of 0. 1 minutes/epoch. The near equivalence between these two numbers - the 
1 -sigma constraint from the observed durations with assumed-flat transit duration variation and 
the numerically estimated rate assuming the 1 -sigma value of mutual inclinations - demonstrates 
our claim. 

3 Bulk composition of the Kepler-36 planets 

3.1 Bayesian Constraints on Kepler-36b's Bulk Composition 

In this section we elaborate briefly on the analysis constraining Kepler-36b's composition. 

Following the Bayesian approach outlined in (83), we account for the (correlated) observa- 
tional uncertainties in the planet mass and radius in a rigorous way when constraining the planet 
composition. We adopt a uniform prior on the planet mass and a flat prior on planet composi- 
tion. The likelihood function (which quantifies how well a specific combination of planet mass 
and composition fit the observed data) is set by the joint M p -R p posterior derived from MCMC 
fits to the transit timing variations and transit light curves (illustrated in Figure 3 in the main 
text). We calculate planet radii at a given mass and composition using the interior structure 
model of (84, 85). 

We first explore a rocky planet scenario for Kepler-36b, in which the planet consists of a 
metal iron core surrounded by a Mgo.gFeo.iSiOa silicate mantle. Our main conclusion is that 
Kepler-36b could be a rocky planet with an Earth-like composition (with ~ 30% of its mass 
in an iron core, and ~ 70% of its mass in a silicate mantle). Marginalizing over planet mass, 
the iron core mass fraction of Kepler-36b is constrained. We find median and 68% credible 
intervals of M corc /M p = 0.29±g;iJ. 

We turn to the possibility that Kepler-36b harbors a water envelope that contributes to its 
transit radius. We follow a similar Bayesian analysis approach as described above, but now 
consider a three component model for the planet structure (iron core, silicate mantle, water en- 
velope) with a flat prior on the fraction of the planet mass in each layer. In this three component 
composition model, inherent degeneracies come to bare on the composition constraints; many 
compositions agree with a single planet mass and radius. The mass of a water envelope on 
Kepler-36b is less than 23% of the total planet mass, at 68% Bayesian confidence. The most 
water rich composition within the 68% credible region that have iron abundances below the (86) 
constraints from silicate collisional stripping simulations is 13% H 2 0. 

Finally, considering the possibility for a H/He gas layer on Kepler-36b, we find that Kepler- 
36b could have no more than 1% of its mass in solar composition H/He layer. 
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Figure S8: Two-parameter joint posterior distributions of the primary model parameters (ex- 
cluding contamination parameters; see §2.2.31) . The densities are plotted logarithmically in 
order to elucidate the nature of the parameter correlations . The indices listed along the diagonal 
indicate which parameter is associated with the corresponding row and column. The parameter 
name corresponding to a given index is indicated in Table IS6l in the "Index" column. 




Figure S9: Posterior distributions in the eccentricity and argument of pericenter planes. 
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0.8 




3D inclination angle (deg) 
Figure S10: Posterior distribution of the mutual inclination, /, of the planetary orbits. 
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8 




Mutual inclination (deg) 

Figure S 1 1 : The change in the transit duration (secular component) over the span of observa- 
tions as a function of mutual inclination for a random sample of 100 initial conditions from the 
posterior. 
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3.2 Bayesian constraints on Kepler-36c's composition 

Kepler-36c has a low enough mean density that it must have a hydrogen-rich envelope (a wa- 
ter envelope is not puffy enough). In the main text (and elaborated here in the supplementary 
online material) we have presented evolution calculations to constrain the H/He envelope mass. 
Assuming a hot-start formation process and passive cooling at the planet's current orbital sepa- 
ration, tight constraints on Kepler-36c's interior entropy and H/He mass obtain. The formation 
mechanism and the full dynamical and insolation histories of the Kepler-36 planets are not 
known, however. 

In this section we complement the planet cooling calculations with a more expansive inves- 
tigation of parameter space, to illustrate the dependence of the planet compositional inference 
on the intrinsic luminosity of the planet. In cases where the planet history is more complicated 
than the isolated cooling calculations have assumed, a systematic offset in the prediction of 
L p (t) could result. The Bayesian contours in Figures IST21 and IS13l offer Kepler-36c composi- 
tion constraints as a function of the planet intrinsic luminosity. 

We consider two scenarios for the composition of the planet's heavy-element interior: a 
rocky (iron and silicate) interior (Figure IS 1 2t . and an ice- rock interior consisting of 60% H 2 
and 40% rock by mass (Figure IS 13b . In both scenarios we consider a range of iron to sili- 
cate mass ratios, and assume a uniform prior on the mass fraction within each compositional 
layer in the planet. Four different choices of Kepler-36c's intrinsic luminosity are shown 
for illustration: L p /M p = 1CT 13 Wkg -1 , a very low intrinsic luminosity to set strong up- 
per limits on the mass of H/He; L p /M p = 10~ 115 Wkg -1 ps L radio /M p , the expected lu- 
minosity from the decay of radiogenic isotopes at 6.8 Gyr assuming bulk Earth abundances; 
L p /M p = 10- 10 - 5 Wkg" 1 « 10L radio /M p ; and L p /M p = 1(T 9 - 5 Wkg^ 1 « 100L radio /M p . 
For low-mass planets with gas layers, the intrinsic luminosity of the planet can have an impor- 
tant effect on the "puffiness" of the gas layer. The H/He envelope mass fraction inferred for 
Kepler-36c is smaller for higher entropy, high L p cases. The hot-start evolution calculations 
favor the higher luminosity scenarios. 

The Bayesian composition constraint contours in Figures IS 121 and IS 1 31 provide a lower 
bound on the uncertainties in Kepler-36c's interior composition and H/He envelope mass at a 
single snapshot of planet intrinsic luminosity. Marginalizing over the planet cooling history to 
account for the uncertainty in L p (t) will further broaden these composition constraints. Uncer- 
tainties quoted in the cooling calculation results are not credible intervals like those presented 
in Figures IST3l and ISTH but are instead estimated by varying one parameter (M p , R p , C p , age, 
incident flux) at a time and adding in quadrature. The heavier computational burden of the 
evolution calculation currently precludes a full Bayesian analysis of the composition uncertain- 
ties. Ultimately, future work will incorporate cooling calculations in the Bayesian approach to 
determine priors on L p (i) and the planet interior entropy. 
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Figure S12: Bayesian constraints on the mass of a solar composition envelope surrounding 
Kepler-36c assuming a rocky interior. The contours plotted are lines of constant probabil- 
ity labeled with a Bayesian confidence value indicating the degree of belief (between and 
1) that the true composition of the planet falls within the contour (given the assumed pri- 
ors). A rocky heavy element interior comprised of an metal iron core and silicate mantle 
is assumed. The Bayesian contours in Figures IS 1 21 and IS 1 3 1 are not marginalized over the 
planet evolution history, but instead show composition constraints for four different choices 
of Kepler-36c's intrinsic luminosity. From high to low H/He content (low to high L p ) the 
values shown are: L p /M p = 1CT 13 Wkg -1 , a very low intrinsic luminosity to set strong up- 
per limits on the mass of H/He; L p /M p = 10" 115 Wkg -1 « L radio /M p , the expected lu- 
minosity from the decay of radiogenic isotopes at 6.8 Gyr assuming bulk Earth abundances; 
L p /M p = 10- 10 - 5 Wkg" 1 « 10L radio /M p ; and L p /M p = 1(T 9 - 5 Wkg" 1 « 100L radio /M p . 
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Figure S13: Bayesian constraints on the mass of a solar composition envelope surrounding 
Kepler-36c assuming an ice-rock interior. This figure is identical to Figure IS 12L but assumes 
the heavy element interior of Kepler-36c is 60% H 2 and 40% rock by mass. The iron to silicate 
ratio abundance is still allowed to vary. 
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3.3 Constraints on Kepler-36c considering thermal evolution 

Kepler-36c is low density, substantially less dense than pure water. As a result, it must have a 
substantial Hydrogen & Helium (H/He) envelope. In order to estimate the planet's mass fraction 
in H/He, we use a coupled thermal evolution and mass loss code (87). These models are based 
on the model used in (88), which has been adapted for modeling Super-Earths following the 
methods of (89). By fully modeling the thermal evolution of the planet's interior, we are able 
to get a much more precise estimate of the the planet's interior structure. Models that only 
compute an instantaneous structure are forced to vary the intrinsic luminosity of the planet over 
several orders of magnitude, which introduces large errors in the composition. 

In particular, the model assumes an adiabatic H/He interior structure that initially has a 
large entropy from formation. We then track the planet's radius up to the present day as it cools 
and contracts. In general, the planet's radius is insensitive to the initial entropy choice by 100 
Myr. The bottom of the H/He adiabat then attaches to an isothermal core with a Earth-like 
rock/iron ratio. Meanwhile, the radiative top of the atmosphere also becomes isothermal once 
this intrinsic flux from this adiabat equals the incident flux at the top of the atmosphere. We find 
cooling rates at the top of the atmosphere, by interpolating between models of the intrinsic flux 
computed on a grid of gravity, Teff, and entropy. These cooling rates were computed using the 
self-consistent, non-gray radiative transfer models described in (29), assuming cloud-free 50x 
solar opacities. 

For the rocky core we use the ANEOS olivine equation of state (EOS) (90) for the rock and 
SESAME FeEOS (91) for the iron. For the H/He envelope we use SCVH (92) and for H 2 we 
use H20-REOS (93). 

Finally, we model the cooling of the rocky core. As the H/He adiabat cools and contracts, the 
isothermal core needs to cool as well. However, because this rock/iron core has a nonzero heat 
capacity, the presence of a large core will slow down the cooling of the atmosphere. Although 
this effect is negligible for gas giants, it becomes dominant for Super-Earths. Models that 
neglect the cooling of the core will predict interiors that are too cold and over predict amount 
of volatiles needed to match the observed radius. We vary the heat capacity of the rock from 
0.5 — 1.0 J/K/g as in (89). We also included radiogenic heating in the rocky core, assuming 
earth like abundances of U, Th, and K. 

In order to determine the current composition of Kepler-36c, we ran thermal evolution mod- 
els without mass loss and adjusted the mass of the core until we matched the observed radius 
at the current age. We then determined the uncertainty in the current composition due to the 
observed uncertainties in the planet's mass, radius, age, and incident flux along with the theo- 
retical uncertainties in the atmospheric albedo, the heat capacity and the rock/iron ratio of the 
core. We varied the albedo from — 0.8 and the iron fraction from pure rock to the maximum 
possible iron fraction from collisional stripping (20). The dominant sources of error were the 
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planet radius, the iron fraction, and the heat capacity of the core, all others were negligible. 
Finally, we were able to get a estimate of the current composition of Kepler-36c of 8.6 ±}; 4 % 
H/He, assuming a water-free interior. 

We also explored the possibility of a water-rich interior for Kepler-36c. To do this we 
inserted a water layer equal in mass to the rocky core in between the rock/iron core and the 
H/He envelope. From the thermal evolution, we find that this water layer should be in the 
molecular and ionic fluid phases and not in high pressure ices. We find that for a water-rich, or 
"Neptune-like" composition, Kepler-36c needs 1.6 ± 0.4% of its mass in H/He. 

Despite the large density contrast seen today, it is possible both planets could have been 
volatile rich in the past. Due their high incident fluxes, both planets are vulnerable to XUV 
driven atmospheric escape of light gases (94, 95). Although it is likely rocky today, we find that 
Kepler-36b could have been as much as 40% H/He when the system was 100 Myr old and the 
stellar XUV flux was ~ 200 times higher (96). With similar models, we find that Kepler-36c 
would have been 35 ± 9% H/He, assuming a dry interior. 

4 Short-term behavior and secular characteristics 
of Kepler-36 

4.1 Evolution of the best-fit model 

The orbital elements show oscillatory patterns as a function of time. Figure IS 141 shows the 
eccentricity vectors plotted versus time. The expected variation in eccentricity is ~ 0.01, which 
is indeed the case for the best-fit model. The eccentricity variation of the smaller planet is larger 
by a factor of the mass ratio, ps 1.8, which is expected due to conservation of momentum during 
the conjunctions. The orbital elements stay approximately constant in between conjunctions, 
then show a jump as they pass through conjunction. As the orientations of the conjunctions 
circulate with time, the orbital elements oscillate. When the components of the eccentricity 
vectors are plotted versus one another, Figure IS 151 the pattern becomes rather intricate. The 
orbital elements remain nearly constant between conjunctions creating the knots in each petal. 
During conjunction the planets undergo the largest change in orbital elements, forming the 
wiggles connecting the knots. The conjunctions drift in inertial space over « 5 conjunctions; 
however, the planets have slightly different orbital elements after 5 conjunctions, so the petals 
drift and eventually form the entire "flower." 

The planet orbits precess with time; Figure IST61 shows the inclination angle and sky nodal 
angle for one of the best-fit models. The timescale of precession in this case is ~ 2.5 x 10 4 
days. 
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Figure S14: Plot of eccentricity vectors for the two planets as a function of time. Kepler-36b is 
plotted in blue, while Kepler-36c is plotted in red. 
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Figure S15: Plot of the two components of the eccentricity vector versus one another. Kepler- 
36b is plotted in blue, while Kepler-36c is plotted in red. 
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Figure S 16: Plot of inclination angles for the two planets and nodal angles as a function of time. 
Kepler-36b is plotted in blue, while Kepler-36c is plotted in red. 
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Figure S17: Predicted transit timing variations based on the models drawn from the posterior 
distribution. The range indicates the 68% confidence range predicted for planet b (blue) and 
planet c (red). 

4.2 Predicted Ephemerides and transit parameters 

Tables IS8l and IS9lprovide predicted times of transit, impact parameters and duration for the next 
year (starting near May 1, 2012) with errors. Figure IS 1 8 I plots predictions over a much longer 
time range after subtracting a linear model. Figure IS 1 8 1 shows the variations in the transit 
durations during that same time. 

5 Long-term stability of solutions 

To assess the long term stability of the system, we selected 10,001 random draws from the pos- 
terior distribution of parameters (masses, initial positions, and initial velocities) and integrated 
each of these for 6.845 x 10 5 years (over 15 million orbits of planet c). We used a symplectic n- 
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Epoch, E 


T transit -((To) + £(P))(hr) 


Impact Parameter 


Duration (hr) 





0.98 ±0.22 


0.130 ±0.078 


6.783±0.077 


1 


0.10 ±0.22 


0.130 ±0.078 


6.783±0.077 


2 


-0.32 ±0.22 


0.130 ±0.078 


6.809±0.077 


3 


-0.47 ±0.21 


0.130 ±0.078 


6.807±0.077 


4 


-0.67 ±0.21 


0.130 ±0.078 


6.808±0.077 


5 


-0.78 ±0.21 


0.130 ±0.078 


6.808±0.077 


6 


-0.87 ±0.21 


0.130 ±0.078 


6.808±0.077 


7 


-1.01 ±0.21 


0.130 ±0.078 


6.808±0.077 


8 


-1.23 ±0.21 


0.130 ±0.078 


6.807±0.077 


9 


-0.10 ±0.19 


0.129 ±0.079 


6.856±0.078 


10 


0.02 ±0.19 


0.129 ±0.079 


6.855±0.078 


11 


0.15 ±0.19 


0.129 ±0.079 


6.856±0.078 


12 


0.37 ±0.19 


0.129 ±0.079 


6.856±0.078 


13 


0.59 ± 0.20 


0.129 ±0.079 


6.856±0.078 


14 


0.75 ±0.21 


0.129 ±0.079 


6.855±0.078 


15 


0.85 ±0.22 


0.129 ±0.079 


6.855±0.078 


16 


2.31 ±0.23 


0.129 ±0.081 


6.872±0.079 


17 


1.78 ±0.25 


0.129 ±0.081 


6.872±0.079 


18 


1.28 ±0.26 


0.129 ±0.081 


6.873±0.079 


19 


0.86 ±0.28 


0.129 ±0.081 


6.873±0.079 


20 


0.44 ±0.30 


0.129 ±0.081 


6.873±0.079 


21 


-0.06 ±0.32 


0.129 ±0.081 


6.872±0.079 


22 


-0.59 ±0.34 


0.129 ±0.081 


6.873±0.079 


23 


-0.51 ±0.37 


0.128 ±0.081 


6.853±0.078 


24 


-1.48 ±0.38 


0.128 ±0.081 


6.853±0.078 


25 


-2.38 ±0.39 


0.128 ±0.081 


6.853±0.078 



Table S8: Predicted ephemerides, impact parameters and durations for Kepler-36b for one year. 
In the table, times are given relative to a linear ephemeris with (T ) = 2456049.39050 (BJD) 
and (P) = 13.85940 days. Epoch E = is near May 1, 2011. 
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Epoch, E 


T transit -((T ) + £(P))(hr) 


Impact Parameter 


Duration (hr) 





-0.81 ±0.05 


0.096 ±0.062 


7.438±0.019 


1 


-0.18 ±0.05 


0.096 ±0.062 


7.438±0.019 


2 


0.19 ±0.04 


0.095 ±0.062 


7.423±0.020 


3 


0.33 ±0.04 


0.095 ±0.062 


7.424±0.021 


4 


0.46 ± 0.03 


0.095 ±0.062 


7.423±0.021 


5 


0.53 ±0.03 


0.095 ±0.062 


7.423±0.021 


6 


0.62 ±0.03 


0.095 ±0.062 


7.424±0.021 


7 


0.77 ± 0.03 


0.095 ±0.062 


7.424±0.021 


8 


0.11 ±0.03 


0.093 ±0.063 


7.395±0.024 


9 


0.03 ±0.04 


0.093 ±0.063 


7.395±0.024 


10 


-0.07 ±0.05 


0.093 ±0.063 


7.394±0.024 


11 


-0.23 ±0.06 


0.093 ±0.063 


7.394±0.024 


12 


-0.36 ±0.07 


0.093 ±0.063 


7.395±0.024 


13 


-0.42 ± 0.09 


0.093 ±0.063 


7.395±0.024 


14 


-1.32 ±0.10 


0.092 ±0.064 


7.384±0.027 


15 


-0.93 ±0.11 


0.092 ±0.064 


7.384±0.027 


16 


-0.59 ±0.12 


0.092 ±0.064 


7.384±0.027 


17 


-0.28 ±0.14 


0.092 ±0.064 


7.384±0.027 


18 


0.06 ±0.15 


0.092 ±0.064 


7.384±0.027 


19 


0.45 ±0.17 


0.092 ±0.064 


7.384±0.027 


20 


0.47 ±0.18 


0.091 ±0.064 


7.396±0.028 


21 


1.17 ±0.19 


0.091 ±0.064 


7.396±0.028 



Table S9: Predicted ephemerides, impact parameters and durations for Kepler-36c for one year. 
In the table, times are given relative to a linear ephemeris with (T ) = 2456044.91461 (BJD) 
and (P) = 16.22407 days. Epoch E = is near April 27, 2011. 



55 



t — i — i — i — i — i — i — i — i — | — i — i — i — i — i — i — i — i — i — | — i — i — i — i — i — i — i — i — i — | — i — i — i — i — i — i — i — i — r 




i i i i i i i i i I i i i i i i i i i I i i i i i i i i i I i i i i i i i i i 

MO 4 2«10 4 3-1 4 4»10 4 

Time (d) 

Figure SI 8: Predicted transit duration variations based on the models drawn from the posterior 
distribution. The range indicates the 68% confidence range predicted for planet b (blue) and 
planet c (red). 
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body map (97, 98) in Jacobi coordinates, with the third-order symplectic corrector developed by 
J. Chambers and reported by Wisdom (2006). In separate work we determined that for the mass 
ratios of Kepler-36b and c (relative to the mass of the star), the Chambers corrector provides the 
optimal balance between precision and computational cost. 

The integrations used a time step of 0.125 days. With the symplectic corrector, the relative 
energy error of the integrations was typically 1CT 12 or smaller, with a negligible secular trend. 
With this time step size, the conjunctions are resolved well. Furthermore, given the degree of 
conservation of the total energy, we believe that this time step is small enough to integrate the 
orbits accurately, despite the relatively close encounters that occur. To test this, we integrated a 
small sample of the initial conditions with a Bulirsch-Stoer integrator and compared the result 
with that from the symplectic n-body map. We confirmed that the output of the two methods 
agrees for short integrations. 

Out of the full set of 10,001 initial conditions, no initial configuration which satisfied the 
Hill criterion was found to have a close encounter, which we defined as occurring when the 
planets were within one Hill radius (99) of each other. We selected a random sample of 100 
initial conditions and integrated these for ~ 140 Million years. Again, there were no close 
encounters. 

Out of the full set of 10,001 initial conditions, 1,171 satisfied (p/a) — (p/a) crit < 1CT 4 . 
Of this subset, 313 were formally Hill stable and 858 did not satisfy the Hill criterion [(p/a) — 
(p/a) crit < 0, (12)] . In other words, these are the initial conditions that do not, or nearly do 
not, satisfy the Hill stability criterion. We integrated these for 13.69 x 10 6 years, again with a 
time step of 0.125 days and with a close encounter threshold of 1.0 Rh- No initial condition 
that satisfied the Hill stability criterion underwent an encounter during the integrations. This 
both further confirms the criterion and supports the validity of our integrations. 

The majority of the 858 initial conditions that do not satisfy the Hill criterion have close 
encounters: 616 passed within one Hill sphere of each other during the 13.69 million years in- 
tegrations. No initial conditions had a close encounter earlier than 10 6 days into the integration. 
Although the Hill stability criterion makes no formal statement about the stability or instability 
of orbits that do not satisfy the criterion, Gladman (1993) (13) found that systems with small ec- 
centricities that do not satisfy the stability criterion typically undergo close encounters on short 
time scales. Since the data constrain eccentricities to be relatively small, we might expect that 
the remaining initial conditions that do not satisfy Hill stability will undergo encounters during 
longer integrations. There is no pressing need to follow up on these exhaustively, though, as 
the goal of this stability analysis was not to understand completely the stability of every initial 
condition. Instead, we have established that ~ 91% of the posterior distribution is long lived 
and does not undergo close encounters. Furthermore, we do not believe that the presence of 
dynamically unstable orbits in the posterior reflects on the stability of the best fit solution, as 
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the Hill unstable portion of the posterior lies on the fringes of the distribution (see Figure |S9T). 

Although the vast majority of the sample of initial conditions are Hill stable, this stability 
criterion alone does not imply that the orbits are Lagrange stable (i.e. that their semimajor axes 
are bounded). The condition of Lagrange stability can further constrain the posterior, just as the 
condition of Hill stability does. This will be explored in future work. 

6 Comparison with other planet systems 

Kepler-36 is an outlier in both separation and difference in density: it has the smallest fractional 
separation of any two adjacent planets with measured masses and radii, and it has the largest 
density contrast, save for Kepler- 18b,c (Figure IS20l ). 

The planets Kepler-36 appear to be located at a border between super-Earths and mini- 
Neptunes: Kepler-36b is the coolest super-Earth with a measured density known to date, while 
Kepler-36c is the hottest mini-Neptune (IS21I) . In the figure we have indicated a suggested 
division between super-Earths and mini-Neptunes at a density of 3.5 g/cc; planets less than this 
density likely have a significant H/He envelope. Any low-density planets (i.e. mini-Neptunes) at 
small orbital separations should be easily detected as these have a higher transit probability and 
deeper transit depth. It is possible that these planets were formed here, but evaporation removed 
the H/He envelopes, causing these planets to evolve towards super-Earths; we have indicated 
this with an arrow in the figure. However, it is quite possible that the hot mini-Neptune region 
will be filled in with further observations. 

7 Availability of data 

The Kepler light curve data used in this analysis can be downloaded from the Mikulski Archives 
for Space Telescopes at the address below by searching for 1 1401755 in the "Kepler ID" field: 

\protect \vrule widthOpt\protect\href{http: //archive. stsci. edu/kepler/ data. 

The spectra used in this analysis have been provided (as ASCII text files) as an attachment 
to this supplement; 

McDonald co-added spectrum : 

1223269s2 . txt 

HIRES spectrum: 

1223269s3 . txt 
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Figure S19: Hill and numerical stability diagram for the 10 4 sets of initial conditions (ICs) 
drawn from the posterior distribution of the photo-dynamical model. The vertical axis is the 
dimensionless analytic Hill- stability criterion due to Marchal & Bozis (1982) (72): positive 
values are Hill stable; horizontal axis is the index for the ICs we tested. The black dots indicate 
ICs that were stable after numerical integration for 2 10 steps of 0.125 d per step (6.84 Myr). The 
colored dots indicate the ICs that became unstable via close approach: blue became unstable 
after > 3 x 10 s time steps, yellow went unstable between 3 x 10 7 — 3 x 10 8 time steps, and red 
went unstable in < 3 x 10 7 time steps. 
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Figure S20: Distribution of fractional separation difference (a 2 /ai — 1) and density ratios 
(pil P\) for adjacent planets in our Solar System (blue) and in Kepler systems (black). Error 
bars are omitted since they are smaller than the size of the point for the Solar System, and are 
difficult to estimate for the Kepler systems (with the exception of Kepler-36). The labels use 
the first letter of the solar system planets, and the Kepler planet number, (i.e. 1 le,b plots the 
ratio of Kepler- 1 If to Kepler 1 le). Kepler-36 is plotted as a red point. 
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Figure S21: Distribution of densities and effective temperatures for super-Earths and mini- 
Neptunes (M p < 10M e ). We have indicated the lack of hot mini-Neptunes (p < 3.5 g/cc, 
T < 1200K) with a dashed blue box. 
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We have attached the 10,001 initial conditions used in the stability analysis (£0) (as an ASCII 
file): 

1223269s4 . txt 

We have attached a draw of 10,000 model parameters from our full posterior (as an ASCII 
text file - with FITS encoding for easier input): 

1223269s5 . txt 

The authors welcome requests for additional information regarding the material presented 
in this paper. 
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